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ABSTRACT 

Miniaturized optical resonators with spatial dimensions of the order of the wavelength of the trapped light 
offer prospects for a variety of new applications like quantum processing or construction of meta-materials. 
Light propagation in these structures is modelled by Maxwell's equations. For a deeper numerical analysis one 
may compute the scattered field when the structure is illuminated or one may compute the resonances of the 
structure. We therefore address in this paper the electromagnetic scattering problem as well as the computation 
of resonances in an open system. For the simulation efficient and reliable numerical methods are required which 
cope with the infinite domain. We use transparent boundary conditions based on the Perfectly Matched Layer 
Method (PML) combined with a novel adaptive strategy to determine optimal discretization parameters like the 
thickness of the sponge layer or the mesh width. Further a novel iterative solver for time-harmonic Maxwell's 
equations is presented. 
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1. INTRODUCTION 

With the advances in nanostructure physics it has become possible to construct light resonators on a lengthscale 
equal to or even smaller than optical wavelengths. 1 ' 2 These nanostructures are large on the atomic scale, 
therefore they can be of complex geometry and they may possess properties not occuring in nature, like an 
effective negative index of refraction 3 which allows in principle to overcome limits in the resolution of optical 
imaging systems. 4 

The numerical simulation of light fields in such structures is a field of ongoing research. In this paper we 
report on finite element methods for the efficient computation of resonances and light propagation in arbitrarily 
shaped structures embedded in simply structured, infinite domains. Section 2 introduces our concept of discretiz- 
ing exterior infinite domains. Section 3 recapitulates a formulation of Maxwell's equations for time-harmonic 
scattering problems. Section 4 introduces an adaptive method for the efficient discretization of the exterior 
domain based on the PML method introduced by Berenger. 5 Section 5 shows the weak formulation of Maxwell's 
equations which is needed for the finite-element method. In Section 6 we shortly introduce a new preconditioner 
for the numerical solution of indefinite time-harmonic Maxwell's equations. Finally, in Sections 8 and 9 we test 
our algorithms on nano-optical real world problems: the computation of resonances and scattering in arrays of 
split-ring resonators and in isolated pyramidal nano-resonators. 
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Figure 1. Infinite domain. The interior domain (left) may contain nearby arbitrarily shaped objects. The exterior 
domain consists of prisms attached to triangular boundary faces of the interior domain. Each prism is the image of the 
unit prism (right) under a bilinear mapping such that the triangles with £ = const are mapped to parallel triangles. For 
each infinite prism we assume constant material parameters. 



In this section we explain how to specify an infinite geometry such that it fits well to the finite element method 
(FEM). The geometry is split into a bounded interior domain f2i nt and an unbounded exterior domain f2 cx t. The 
interior domain may contain nearby arbitrary shaped structures such as spheres or thin layers. The geometry 
in the exterior domain is more restricted. However, the construction we propose is general enough to deal with 
typical geometries of optical devices. 

We assume that the boundary T of the interior domain consists of triangles. A boundary triangle FcT 
is called transparent if F C <9f2 xt- Further we introduce the unit prism P n = {(771,772,0 € R 3 : V1tV2,£, > 
0,771 + r)2 < 1}. An exterior domain is admissible if the following conditions are satisfied. For each boundary 
triangle F there exists a bilinear one-to-one mapping Qp from the unit prism into the exterior domain f2 cx t such 
that each triangle T p = 172, £) C P u : £ = p} is mapped onto a triangle parallel to the face and such that 
the bottom triangle of the unit patch is mapped onto the corresponding face, QfTq = F, cf. Figure 1. The 
image of Qp is denoted by Pp. Hence we attach the infinite prism Pp to the transparent face F. It must hold 
true that Jl ex t = Up Pp. Further we demand the following matching condition. If Qf(t?i, 772, 0) = Qf> , 772 1 0) , 
that is F and F' have a common point, then Qf(?7i,?72,£) = Qf'(t1i> W^i f° r au £ ^ R+- 

The surface S p = UpQpT p looks like a "stretched" transparent boundary of the interior domain. Hence the 
coordinate £ is chosen consistently for all prism such that it serves as a generalized distance variable. This is 
essential for the pole condition concept developed by Frank Schmidt. 6 For later purposes we introduce the 
truncated unit prism P p = {(?7;l, »72; £) € R 3 : 771, 772 > 0, 771 + 772 < 1, < £ < p} and the truncated exterior 
domain Q p = \JpQpP p . 

If there exist triangles on T which are not transparent, then either boundary conditions must be imposed on 
them, or they must be identified with other periodic triangles (e.g., when f2 is a cell of a periodic structure). 
However for simplicity we assume that S7i nt U fl ox t — R 3 in rest of the paper. 



3. SCATTERING PROBLEMS 

Monochromatic light propagation in an optical material is modelled by the time-harmonic Maxwell's equations 



2. GEOMETRIC CONFIGURATION 



curl fi 1 (x) curl E (x) - u?e (x) E (x) = 0, 
div e (x) E (x) = 0, 



(la) 
(lb) 



which may be derived from Maxwell's equations when assuming a time dependency of the electric field as 
E(x, i) = E (x) exp(— iuut) with angular frequency w. The dielectric tensor e and the permeability tensor fi 
are L°° functions of the spatial variable x = (xi, X2, X3). In addition we assume that the tensors e and /j arc 
constant on each infinite prism as defined in the previous section. For simplicity assume that the dielectric and 
the permeability tensors are isotropic so they may be treated as scalar valued functions. Recall that any solution 
to (la) with lo ^ also meets the divergence condition (lb). 

A scattering problem may be defined as follows: Given an incoming electric field Ei nc satisfying the time- 
harmonic Maxwell's equations (1) for a fixed angular frequency oj in the exterior domain, compute the total 
electric field E satisfying (1) in £l; n t U fi C xt, such that the scattered field E sc = E — Ei nc defined on £l ext is 
outward radiating. For a precise definition of when a field is outward radiating we refer to Schmidt. 6 Hence the 
scattering problem splits into an interior subproblem for Ei nt = E|Q. nt on f2; nt 

curl ^i _1 curl E int — w 2 eEi nt = 0, 

and an exterior subproblem on f2 ext 

curl /j^curl E sc — a> 2 eE sc = 0. 

These subproblcms are coupled by the following matching conditions: 

Eint x n = (E inc + E sc ) x ft 
/j^curl Ei n t x n = /i _1 curl (E inc + E sc ) x n 

on the boundary d^mt- 

4. ADAPTIVE PML METHOD 

The perfectly matched layer method was originally introduced by Berenger in 1994. 5 The idea is to discretize 
a complex continued field in the exterior domain which decays exponentially fast with growing distance to the 
interior-exterior domain coupling boundary. This way a truncation of the exterior domain only results in small 
artificial reflections. The exponential convergence of the method with growing thickness of the sponge layer was 
proven for homogeneous exterior domains by Lassas and Somersalo. 7 ' 8 An alternative proof with a generalization 
to a certain type of inhomogeneous exterior domain is given by Hohage et al. Nevertheless as shown in our 
paper 10 the PML method intrinsically fails for certain types of exterior domains such as layered media. This 
is due to a possible total reflection at material interfaces. In this case there exists a critical angle of incidence 
for which the resulting field in the exterior domain is neither propagating nor evanescent. Here we show that 
it is possible to overcome these difficulties when using an adaptive method for the discretization of the exterior 
domain problem. We assume the following expansion of the scattered field in the exterior domain 

Esc(vuV2,0~ [ c( m , m ,a)e ik ^da (4) 

with ^Rk^(a) > 0, ^sk^(a) > and a bounded function 0(771,772,0;). Hence E sc is a superposition of outgoing or 
evanescent waves in £ direction. In our notation we have assumed that there exists a global (771, 772, ^-coordinate 
system for the exterior domain. But in the following only the global meaning of the £ coordinate as explained 
in Section 2 will be used, so 771 and 772 may also be considered as coordinates of a local chart for a subdomain of 
dQint- For 7 = 1 + ia the complex continuation, £ 7^, E SCj7 (-, •,£;) = E sc (-, -,7£) gives 

|E SC , 7 | < e- KX *C, (5) 

with k = min a {3fc(a), adik(a)}. Therefore E SCi7 (-,-,£) decays exponentially fast with growing generalized 
distance £ to the coupling boundary. The idea is to restrict the complex continuation of the exterior domain 
problem to a truncated domain fi p and to impose a zero Neumann boundary condition at dVl p . In the next 
section we will give a corresponding variational problem which can be discretized with the finite element method 
where we will use a tensor product ansatz in the truncated exterior domain il p based on the triangulation of the 



(2) 
(3) 



Algorithm 1 Adaptive PML method 



Require: e, a, h int , K min 

Compute iV p . w and £ max depending on h- lnt and finite element order 
while (not converged) do 
£ =0.0;<£i = /i int ;JV=l; 
while (-ln(e)/(^jvcr) < K min ) do 

S.n+1 = + max{/i int , 2tt<t£ n / (- ln(e))/7V p . w }. 
if (Cat+i > 1/e) then 

break 
else 

N = N+ 1 
end if 
end while 

Compute solution u with PML discretization {£o,£i> ■ • • ,£n} 
if < then 

converged 
else if £ w > £max then 

break 
else 

^min 

end if 
end while 



surface <9fii n t and a ID mesh in ^-direction, {0, £1,^2, • ■ • £jv}- hi this section we present an algorithm for the 
automatic determination of optimal discretization points £j . 

As can be seen from Equation (5) the PML method only effects the outgoing part with Sftfc^ strictly larger 
than zero. Field contributions with an large 3?fc^ component arc efficiently damped out. Furthermore evanescent 
field contributions are damped out independently of the complex continuation. For a proper approximation of 
the oscillatory and exponential behavior a discretization that is fine enough is needed to resolve the field. In 
contrast to that anomalous modes or "near anomalous" modes with ~ enforce the usage of a large p but 
can be well approximated with a relatively coarse discretization in £. Such "near anomalous" modes typically 
occur in the periodic setting but may also be present for isolated structures with a layered exterior domain. 11 
Hence for an efficient numerical approximation of the scattered field one must use an adaptive discretization. 
It is useful to think of the complex continuation as a high-frequency filter. With a growing distance £ to the 
interior coupling boundary the higher frequency contributions are damped out so that the discretization can be 
coarsened. 

For a given threshold e we introduce the cut-off function 

Kco, e (O = -M«0/£ • 

At £' > each component in the expansion (4) with k > K C o,e(£') is damped out by a factor smaller than the 
threshold e, 

e -<' < e -«co,«(«')€ = e Me) = ^ 

Assuming that this damping is sufficient we are allowed to select a discretization which must only approximate 
the lower frequency parts with k < K co ,e(£) for £ > If we use a fixed number iV p . w of discretization points 
per (generalized) wavelength 2tt/k we get the following formula for the a priori determination of the local mesh 
width /i(£) = 27r<7//t C o,e(£)/-/Vp. w . Since K co ,e(£) — > 00 for £ — > the local mesh width is zero at £ = 0. As it is not 
reasonable to use a finer discretization in the exterior domain than in the interior domain we bound the local 
mesh width by the minimum mesh width h m t of the interior domain discretization on the coupling boundary 

h(£) = max{/i int , 27rcr/K CO!e (£)/7Vp. w }. 



■'out 



Xl 



Figure 2. Test problem for adaptive PML discretization. A plane wave is incident under an angle i? from the lower 
material with refractive index n su b = 1.5. The upper material consists of air (n sup = 1.0). According to Snell's law the 
field is totally reflected for an incident angle greater or equal to the critical angle # c = 180 ■ asin(1.0/1.5)/7r ~ 41.81. 
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Figure 3. Left: Field energy error in the interior domain. The three data sets (o, A, +) correspond to different refinement 
levels of the interior domain. Right: Zoom into the left figure near the critical angle. 



The parameters e and iV p w are also fixed accordingly to the interior domain discretization quality. The grid 
{£o; £ij £a> • • • } is recursively constructed by 

This way £ n grows exponentially with n. To truncate the grid we assume that components in the expansion with 
k < Kmin can be neglected so that the grid {£o,£i, ■ ■ ■ ,?at} is determined by k C o,<;(£,n) < K m i n < n C o,e(£,N-i)- 
As an a posteriori control we check if the field is indeed sufficiently damped out at £/v, ||u(-,£jv)|| < 
Otherwise we recompute the solution with K m i n — > K m i n /2 *. Since for an anomalous mode the field is not damped 
at all we restrict the maximum £/y to £jv < ft/ko/e. The pseudocode to the algorithm is given in Algorithm 1. 

To demonstrate the performance of the adaptive PML algorithm we compute the reflection of a plane wave 
at a material interface, cf. Figure 2. In x%— direction we use Bloch periodic boundary conditions. 12 We vary 
the angle of incidence from = 20° to i9 = 60°. Further the incoming field is rotated along the 23 axis by an 
angle of 45°, so that the incidence is twofold oblique (conical). Hence the unit direction of the incoming field is 
equal to fc = (cos 45° sin$, cos$, sin 45° sin$). We use an interior domain of size 1.5 x 1 in wavelength scales. To 
measure the error we compute the field energy within the interior domain and compare it to the analytic value. 
In Figure 3 the error is plotted for different refinement levels of the interior domain. The "+" line corresponds 
to the finest level. In Figure 4 the automatically adapted thickness of the PML is plotted (left) and the number 
of discretization points N in £ direction (right). As expected a huge layer is constructed automatically at the 
critical angle, whereas the total number of discretization points remains moderate. As can be seen in Figure 3 

"This strategy proved useful in many experiments. However we consider to refine it. 
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Figure 4. Left: Thickness of the PML layer in unit lengths. At the critical angle the thickness is up to 10 4 times larger 
than the size of the interior domain. Right: Number of discretization points £j used in the radial direction (£2). Although 
the required thickness of the layer is huge the number of unknowns used in the PML layer remains moderate. 
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0.359850 


0.335129 
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0.159358 


0.166207 
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0.048779 


0.049502 
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0.012911 


0.012912 
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0.003274 


0.003266 


5 


0.000205 


0.000820 


6 


0.000206 


0.000205 


7 


0.000051 


0.000051 



Table 1. Convergence of field energy at the critical angle of incidence. The first column corresponds to the interior 
mesh refinement step. The relative error of the electric field energy in the interior domain is given in the second column. 
AE — |||E e x||^2 — || Eh || 2 L 2 1 / 1| E ex ||^2- The third column displays the relative error of the magnetic field energy AE' = 
1 1 1 curl E e x||^2 — 1 1 curl E h ||^2 |/||curl E ex || is given. For fixed PML thickness the solution converges towards the analytical 
result as the interior mesh is refined. 



the maximum error appears at the critical angle. From that one may suspect a failure of the automatic PML 
adaption. But a closer analysis reveals that the chosen discretization in the PML layer is sufficient as can be seen 
from Table 1. Here the thickness of the perfectly matched layer has been fixed and we further refined the interior 
domain. By this means we observe convergence to the true solution but the convergence rate is halved at the 
critical angle. Hence the maximum error at the critical angle is caused by an insufficient interior discretization. 
We conjecture that this is due to a dispersion effect. Since near the critical angle the wave E out is traveling 
mainly along the x\— direction it reenters the periodic domain, leading to large "path length". 



5. VARIATIONAL FORMULATION 

So far the overall scattering problem was given as an interior domain problem coupled to an exterior domain 
problem via boundary matching conditions. In this section we give (without proof) a variational problem in 
H (curl , fiUOext) f° r the computation of the composed field E with E = Ej nt in £l; nt and E = E SCi7 + n(E; nc x n) 
in fiext • Details for the 2D case are given in our paper. 13 Here n is the extension operator defined as n(Ei nc xn) = 
X[o,e)(l-£A)(Einc x ft). 

For each face F of the transparent boundary Jp(r]i, 772, £) denotes the Jacobian of the mapping Qf(vii V2,C)- 
Further we introduce the pulled back field 11,(771, 772, £) = J t \i{Qp(rji, 772, £)) for any field defined on fi ex t- With 
the definition 

curl 7 = (d V2 dt,-d m + -<%, d m - d m ) 

7 7 

and the transformed tensors /j* = | J|J _1 /iJ _t and £* = |J|J _1 £j _t the composed field E satisfies 

/ curl */i _1 curl E - w 2 *£E + 7 V] / curl 7 **7J^ 1 curl 7 E» - w 2 **£„E* = 

■/flint F JPu 

— <& 7i _1 curl 3 E inc x 77 + 7^] / curl 7 *,/jJ 1 curl 7 n(E inc x ft) — w 2 \l>*£n(Ei nc x ft). 

Although this equation looks complicated it can be easily discretized with finite elements. In fact the terms 
in the sum over the faces F are already given in unit coordinates of the prism. So it is advantageous to use a 
prismatoidal mesh in the exterior domain. This way we fix the global discretization points {£0 = 0, £1, • • • ,£jv} 
as described in the previous section and split the truncated exterior domain into the prisms Qf{Ps,i +1 \ 
with i < N. In the interior domain we use a tetrahedral mesh which we fit non-overlapping to the exterior domain 
mesh. Introducing the bilinear forms 



flint (*, *) 


= / curl *tj ^urlE 

Jn iat 


knt (*, *) 


= [ *£E 

J flint 


fl 7 (*,*) 


= 7 V / curl 7 *»tj^ 1 curl 7 E» 
F Jp p 




= 7V / **£»E* 
F Jp p 


a(*,*) 


= flint (*,*) + 7 (*, *) 


&(*,*) 


= &int(*,*)+M*,*) 



and 

9 (*) = - / * 77 _1 curl 3 Ei„ c x 77 

• / *>int 

the variational problem truncated to f2; nt U VL p can be casted to 

a (*, E) - u?b (*, E) = g (*) + a 7 (*, n(E inc x rl)) - uj\ (*, n(E inc x ft)) 



for all <G H (curl , f2i nt U ft p ). To discretize this variational problem we use Nedelec's vectorial finite elements 14 
with local ansatz functions {vi, v 2 , . . . ,v„}. Making the ansatz E = ^UiVi this yields the algebraic system 

(A - cj 2 )Bu = f 

with fi = g (vi) + a 7 (v i; n(E inc x n)) - oj 2 b 7 (vi, n(E inc x n)) , A itj = a (v ; , vj) and B iyj accordingly 

6. TIME DOMAIN PRECONDITIONER 

In this section we propose a novel preconditioner for the algebraic system 

(A - cj 2 B)u = f (6) 

derived in the last section. Since u> 2 > this system is indefinite. Hence standard multigrid methods will suffer 
from slow convergence rates or may even not converge. Other numerical methods like the Finite Difference 
Time Domain method do not start from the time-harmonic Maxwell's equations. Instead they simulate temporal 
transient effects. For practical purposes the computation time is prohibitively large until the steady state is 
reached. 15 Even worse, the usage of an cxplicite time stepping scheme forces the usage of very small time steps 
to avoid instabilities. 

Here we propose a preconditioner for the time-harmonic system which makes use of the fact that the solution 
we want to compute is the steady state solution to a transient process. This is the reason why we call this 
preconditioner "time domain" preconditioner. Instead of using time dependent Maxwell's equations one may use 
another dynamical system whose steady state solution is the field in mind. For example the above solution u is 
the steady state solution to the time dependent problem 

i^Bu(t) = Au(t) - fe~ iujH 

which looks like the time dependent Schrodinger equation. But we may also start from a time discrete system, 
such as 

i-(Bu n+1 - Bu n ) = Au n+1 - p n+1 f (7) 

T 

with p = ( i + i^ T ) corresponding to a time discretization of the above Schrodinger like equation with the im- 
plicite Euler method. One proves that if the original system has a steady state solution then the sequence 
{uo,p^ 1 ui,p^ 2 U2 7 ■ ■ ■} converges to the solution u of Equation (6) indepently of the selected time step r. In our 
code we typically fix r = O.I/w 2 . Starting from a randomized initial guess uo we compute a fixed number N of 
iterations to the recursion formula (7) which yields the sequence {mo,P _1 ui, ■ • • ,p~ n un}- In each iteration step 
the arising system is solved by an multigrid method up to a moderate accuracy. We then compute the minimum 
residual solution within the space spanned by the last M < N vectors in this sequence. The so constructed 
approximate solver is used as a preconditioner for a standard iterative method for indefinite problems such as 
GMRES or BCGSTAB. 16 

Other discrete schemes may be used to improve the convergence to the steady state solution. For example 
it is promising to use schemes stemming from higher order Runge-Kutta methods or multi-step methods for the 
discretization of the original wave equation or the Schrodinger like equation above. 

7. RESONANCE PROBLEMS 

A resonance is a purely outgoing field which satisfies the time harmonic Maxwell's equation for a certain ui e C. 
We again assume that an expansion as in Equation (4) is valid but we must drop the assumption ^sk^(a) > 0. 
Hence a resonance mode may exponentially grow with the generalized distance £. In this case we must choose a 
large enough in the PML method to achieve an exponential damping of the complex continued solution. Using 
a finite element discretization as in the previous section we end up with the algebraic eigenvalue problem 



Au = lu 2 Bu. 



(8) 



1 1 — - 
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Figure 5. (a) Visualization of a part of the tetrahedra of the spatial discretization of the interior SRR geometry. Dark gray 
tetrahedra: gold SRR; light gray: air; gray: ITO. Periodic boundary conditions apply at the right/left and front/back. 
Prism elements discretizing the exterior domain (on the top/bottom) are not shown, (b) Transmission spectra of light 
fields incident onto an SRR for different angles of incidence (For details see also 17 ). The transmission minimum at 
A ~ 1.5/xm is due to the excitation of the fundamental resonance of the SRRs. (See original publication for images with 
higher resolution.) 



8. META-MATERIALS: SPLIT RING RESONATORS 

Split-ring resonators (SRR's) can be understood as small LC circuits consisting of an inductance L and a 
capacitance C. The circuit can be driven by applying external electromagnetic fields. Near the resonance 
frequency of the LC— oscillator the induced current can lead to a magnetic field opposing the external magnetic 
field. When the SRR's are small enough and closely packed - such that the system can be described as an 
effective medium - the induced opposing magnetic field corresponds to an effective negative permeability, \i < 0, 
of the medium. 

Arrays of gold SRR's with resonances in the NIR and in the optical regime can be experimentally realized using 
electron-beam lithography, where typical dimensions are one order of magnitude smaller than NIR wavelengths. 
Details on the production can be found in Linden and Enkrich. 1 ' 2 

Due to the small dimensions of the LC circuits their resonances are in the NIR and optical regime. 2 Fig- 
ure 5(a) shows the tetrahedral discretization of the interior domain of the geometry. Figure 5(b) shows results 
from FEM simulations of light scattering off a periodic array of SRR's for different angles of the incident light. 17 
At A <~ 1.5/xm the transmission is strongly reduced due to the excitation of a resonance of the array of SRR's. 
This excitation occurs for all investigated angles of incidence. 

When one is interested to learn about resonances, obviously it is rather indirect and time-consuming to 
calculate the scattering response of some incident light field and then to conclude the properties of the resonance. 
We have therefore also directly computed resonances of SRR's by solving Eqn. (8). Special care has to be taken 
in constructing appropriate PML layers, as in this case the previously described adaptive strategy for the PML is 
more involved. We have therefore set these parameters by hand. We have computed the fundamental resonance 
with u = 1.302023 • 10 15 - 0.399851 • I0 15 i. 




Figure 6. Pyramidal nano resonator. The structure in mounted on a Gallium Arsenid (GaAs) substrate. The middle 
figure shows the field amplitude in the x-z plane. The right figure shows the field amplitude in the y-z plane. (See original 
publication for images with higher resolution.) 

9. PYRAMIDAL NANO-RESONATOR 

This type of nano resonators is proposed to be used as optical element in quantum information processing. 18 
The structure is as depicted in Figure 6 (left). We have simulated the illumination of the structure with a 
plane wave of a vacuum wavelength A = 1.55/xm and a unit direction k = (VoU, 0, — \/0J3). The incident field 
was polarized in x-direction. Figure 6 shows the field amplitude in the computational domain. More than five 
million of unknowns were used in the discretization. With the preoconditioner proposed in Section 6 (N=20) the 
GMRES method exhibited a convergence rate of 0.8. 
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